A Computational Framework for Infinite-Dimensional Bayesian Inverse Problems, Part II: Stochastic Newton MCMC with Application to Ice Sheet Flow Inverse Problems

We address the numerical solution of infinite-dimensional inverse problems in the framework of Bayesian inference. In the Part I companion to this paper, we considered the linearized infinite-dimensional inverse problem in which the mean and covariance of the posterior parameter measure were approximated, respectively, by the maximum a posteriori (MAP) solution and the inverse of the Hessian of the negative log posterior at the MAP. Here in Part II, we relax the linearization assumption and consider the fully nonlinear infinite-dimensional inverse problem. We first introduce the infinite-dimensional setting, and next present a finite-dimensional approximation that is consistent with this infinite-dimensional setting. The resulting high-dimensional finite-dimensional posterior is then explored using a Markov chain Monte Carlo (MCMC) sampling method. Sampling the posterior is challenging due to high dimensionality of the discretized parameter space, and whenever the forward model is governed by expensive-to-solve partial differential equations (PDEs). To address the challenges of sampling high-dimensional pdfs arising from Bayesian inverse problems governed by PDEs, we build on the stochastic Newton MCMC method. This method exploits problem structure by taking as a proposal density a local Gaussian approximation of the posterior pdf that is based on local gradient and Hessian information. Construction and manipulation of the Hessian is made tractable by invoking a low-rank approximation of its data misfit component, motivated by the fact that the data are typically informative about only a low-dimensional manifold in parameter space. However, even with the use of this low-rank approximation, the need to recompute the Hessian at each proposed sample point can be prohibitive for large-scale problems. To address this problem, here we introduce an approximation of the stochastic Newton proposal in which we compute the low-rank-based Hessian at just the MAP point, and then reuse this Hessian at each MCMC step. We compare the performance of the proposed method to the original stochastic Newton MCMC method (with dynamically-computed Hessian) and to an independence sampler that employs as an MCMC proposal a Hessian-based Gaussian centered at the MAP point. These three Hessian-based methods are identical in the case of a linear inverse problem with Gaussian prior and noise models, but differ for a non-Gaussian posterior. The comparison of the three methods is conducted on a synthetic ice sheet inverse problem. The inverse problem is to infer the coefficient field of the boundary condition at the base of the ice, given free surface velocity observations and a nonlinear Stokes model describing the non-Newtonian, viscous, creeping, incompressible flow of the ice. For this problem, the stochastic Newton MCMC method with a MAP-based Hessian converges at least as rapidly as the original stochastic Newton MCMC method (with dynamically-computed Hessian), but is far cheaper since it avoids recomputing the Hessian at each step. On the other hand, it is more expensive per sample than the independence sampler due to its need to compute the gradient at each step; however, its convergence is significantly more rapid, and thus overall it is much cheaper (in terms of forward PDE solves). Therefore, the method proposed in this paper performs the best of the three Hessian-based MCMC methods. Finally, we present extensive analysis and interpretation of the posterior distribution, and classify directions in parameter space based on the extent to which they are informed by the prior or the observations.